---
title: "Power Calculation based on simr simulation"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

#install.packages("simr")
library(simr)
```


```{r}
## building model object from covariates
sub <- factor(1:15) #set number of subjects
condition <- factor(c("baseline","distraction"))
modality <- factor(c("vibrotactile","thermonociceptive","thermonociceptive","vibrotactile"))


sub_full <- rep(sub, 4) #repeat subjects for each condition/modality combination
condition_full <- rep(rep(condition, each=15),2 ) 
modality_full <- rep(modality, each=15)

covars <- data.frame(id=sub_full, condition=condition_full, modality=modality_full) #combine all variables to one dataframe

covars


## Intercept and slopes for modality (from Moulders2020), condition (adapted from Wierda 2010), 
## modality:condition (guessed since no data available)

fixed <- c( 0.296,0.2,-0.7,-0.14)
            
## Random intercepts for participants
rand <- list(0.01)

## residual variance, res <-0.0966923
res <- 0.09349

model <- makeLmer(y ~ modality*condition + (1|id), fixef=fixed, VarCorr=rand, sigma=res, data=covars) #Linear mixed model using the predefined variables
model #to check model
```


```{r}
#create power curve with initial amount of participants. Testing the effect of the interaction, based on sample size and for 100 simulations. Cortex specific sample size guidelines: power needs to reach 90% and alpha=0.02.
p_curve <- powerCurve(model, test=fixed("modality:condition"), along="id", nsim=100, alpha=0.02, seed=123) 
plot(p_curve)
powerSim(model, nsim=100, test = fixed("modality:condition"), alpha=0.02,seed=123) #same power calculation, without plot

#extend model to have a bigger sample size
model_long <- extend(model, along="id", n=25)

#use extended model for Power calculation
p_curve_long <- powerCurve(model_long, test=fixed("modality:condition"), along="id", nsim=25, alpha=0.02,seed=123) 

plot(p_curve_long) #show plot of power calculation result

sim_cond <- powerSim(model_long, nsim=25, test = fixed("modality:condition"), alpha=0.02,seed=123) #power calculation without plot
sim_cond

#powerSim(model, nsim=25, test = fcompare (~condition), alpha=0.02,seed=123) #comparing our model to a simpler model using only the factor condition
# lastResult()$warning #shows the last warning in case of problem with running the package

```


